Review

Bayesian analysis is a contrast to frequentist methods

• The objective of Bayesian analysis is to compute a posterior distribution
  

• Contrasts with frequentist statistics with objective to compute a point estimate and confidence interval from a sample

• Bayesian models allow expressing prior information as a prior distribution

• The posterior distribution is said to quantify our current belief

  • Belief is based on the posterior distribution
    
  • We update beliefs based on additional data or evidence
    
  • A critical difference with frequentist models which must be computed from a complete sample

• Inference can be performed on the posterior distribution by finding the maximum a postiori (MAP) value and a credible interval

• Predictions are made by simulating from the posterior distribution

Review

Two functions must be defined to compute the posterior distribution

• The likelihood for the model being used
  • The likelihood function includes the parameters for the model
    
  • Example; for Binomial likelihood is the probability of success
    
  • Example; for Normal likelihood is the mean, \(\mu\) and variance. \(\sigma^2\)
• The prior distribution encodes the information we have in advance about the model parameters
  • For simple cases is the conjugate distribution
    
  • The posterior distribution is in the conjugate family
    
  • Example, for binomial likelihood, the conjugate is the \(Beta(\alpha,\beta)\)
    
  • Example, for Normal likelihood, the normal distribution is the conjugate for the mean, \(\mu\)

Review

Several solution methods; we apply two possibilities

• Use an analytical solution; posterior in family of conjugate prior

  • Update conjugate prior using evidence (observations or data)

• Grid sample over range of model parameter values

• In both cases, goal is to compute the posterior distribution over the plausible range of parameters

  • Example; for binomial likelihood, range of parameter values, \(0 \le\theta \le 1.0\)
    
  • Example; for Normal likelihood, range of parameter values around empirical values of \(\mu\) and \(\sigma > 0\)

Review

What are the properties of the Beta distribution?

Beta distribution for different parameter values

Review

Consider the product of a Binomial likelihood and a Beta prior

• Define the evidence as \(n\) trials with \(z\) successes
  • These are the actual data or observations for your problem
• Prior is a Beta distribution with parameters \(a\) and \(b\), or the vector \(\theta = (a,b)\)
  • \(a\) is the pseudo count of sussesses
    
  • \(b\) is the pseudo count of failure
    
  • For an uninformative prior, \(U(0,1) = Beta(1,1)\), so \(a=0,\ b=0\), since we have no prior information
    
  • With uninformative prior, posterior distribuytion = likelihood funciton
• From Bayes Theorem the distribution of the posterior:

\[\begin{align} posterior(\theta | z, n) &= \frac{likelihood(z,n | \theta)\ prior(\theta)}{data\ distribution (z,n)} \\ Beta(\theta | z, n, a, b) &= \frac{Binomial(z,n | \theta)\ Beta(a+1, b+1)}{p(z,n)} \\ &= Beta(z + a -1,\ n-z+b-1) \end{align}\]

Review

Consider example with:
- Prior pseudo counts \([1,9]\), successes \(a = 1 + 1\) and failures, \(b = 9 + 1\)
- Evidence, successes \(= 10\) and failures, \(= 30\)
- Posterior is \(Beta(10 + 2 -1,\ 40 - 10 + 10 -1) = Beta(11,\ 39)\)

Prior, likelihood and posterior for distracted driving

Review

What are the 95% credible intervals for \(Beta(11,\ 39)\)?

Ppsterior density of probability of distract drivers

Review

Simulate the probabilities of distracted drivers for the next 10 cars with posterior, \(Beta(11,\ 39)\)?

Probability of distract drivers for next 10 cars

Introduction

How can we extend Bayes models to more complex problems?

• For simple problems we can use a conjugate prior and posterior

• Unlikely posterior distribution will be a simple conjugate

• Need to perform sampling to compute approximation of complex posterior

• We need highly efficient sampling methods for complex problems

Problem with Grid Sampling

Real-world Bayes models have large numbers of parameters, into the millions

• Naive approach is simple grid sampling

  • Sample across dimensions of the parameter space

• Grid sampling does not scale

• Consider this thought experiment, sampling dimension 100 times:

  • 1-parameter model: \(100\) samples
    
  • 2-parameter model: \(100^2 = 10000\) samples
    
  • 3-parameter model: \(100^3 = 10^5\) samples
  • 100-parameter model: \(100^{100} = 10^{102}\) samples

• Computational complexity of grid sampling has exponential scaling with dimensionality

• Need a better approach!

Introduction to Markov Chain Monte Carlo

Large-scale Bayesian models need highly efficient sampling methods

• Markov chain Monte Carlo (MCMC) sampling is efficient and scalable

• Rather than sampling on a grid MCMC methods sample distributions randomly

• Requires effort to understand how the algorithm works and what to do when things go wrong

What is a Markov process?

MCMC sampling uses a chain of Markov sampling processes

• A Markov process is a stochastic process that transitions from a current state, \(x_t\), to some next state, \(x_{t+1}\), with probability \(\Pi\)
  • No dependency on past states
• Summarize properties of a Markov process:
  • Probability of state transition is parameterized by a matrix of probabilities, \(\Pi\), of dim N X N for N possible states
    
  • \(\Pi\) only depends on the current state, \(x_t\)
    
  • Transition can be to current state.

What is a Markov process?

Since Markov transition process depends only on the current state a Markov process is memoryless

• We can express the sequence of a Markov transition processes as:

\[P(X_{t + 1}| X_t = x_t, x_{t-1}, x_{t-2}, \ldots, x_0) = p(X_{t + 1}| x_t)\]

• The Markov process is memoryless
  • Transition probability only depends on the current state, \(x_t\)
    
  • No dependency on any previous states, \(\{x_{t-1}, x_{t-2}, \ldots, x_0 \}\).

What is a Markov process?

For system with \(N\) possible states we can write the transition probability matrix, \(\Pi\):

\[\Pi = \begin{bmatrix} \pi_{1,1} & \pi_{1,2} & \cdots & \pi_{1, N}\\ \pi_{2,1} & \pi_{2,2} & \cdots & \pi_{2,N}\\ \cdots & \cdots & \cdots & \cdots \\ \pi_{N,i} & \pi_{N,2} & \cdots & \pi_{N,N} \end{bmatrix}\\ where\\ \pi_{i,j} = probability\ of\ transition\ from\ state\ i\ to\ state\ j\\ and\\ \pi_{i,i} = probability\ of\ staying\ in\ state\ i\\ further\\ \pi_{i,j} \ne \pi_{j,i}\ in\ general \]

Example of a Markov Process

To make the foregoing more concrete let’s construct a simple example. We will start with a system of 3 states, \(\{ x_1, x_2, x_3 \}\). The transition matrix is:

\[\Pi = \begin{bmatrix} \pi_{1,1} & \pi_{1,2} & \pi_{1,3}\\ \pi_{2,1} & \pi_{2,2} & \pi_{2,3}\\ \pi_{3,1} & \pi_{3,2} & \pi_{3,3} \end{bmatrix} = \begin{bmatrix} 0.5 & 0.0 & 0.6\\ 0.2 & 0.3 & 0.4\\ 0.3 & 0.7 & 0.0 \end{bmatrix} \]

• Some key points to notice
  • The probabilities of transition from a state is given in a column
    
  • The probabilities in each column must add to 1.0
    
  • The probabilities of a transition to the same state are along the diagonal of the matrix
• Some transitions not possible have a probability of 0.0
  • Example, \(\pi_{2,1} = 0\); cannot transition from state 2 to 1
  • Example, \(\pi_{3,3} = 0\); cannot remain in state 3

Example of a Markov Process

Let’s apply a probability matrix to a set of possible states

• The state vector represents being in the first state at time step \(t\); \(\vec{x_t} = [1,0,0]\)
  
• After a state transition, we compute the probability of being in each of the three possible states at the next time step, \(t+1\):

\[\vec{x}_{t+1} = \Pi\ \vec{x}_t = \begin{bmatrix} 0.5 & 0.0 & 0.6\\ 0.2 & 0.3 & 0.4\\ 0.3 & 0.7 & 0.0 \end{bmatrix} \begin{bmatrix} 1\\ 0\\ 0 \end{bmatrix} = \begin{bmatrix} 0.5 \\ 0.2 \\ 0.3 \end{bmatrix} \]

From Markov process to Markov chain

The foregoing is a single step of a Markov process

• What happens when there is a series of transitions?

• A sequence of such transitions is known as a Markov chain

• There are two major behaviors observed with Markov Chains

• Episodic Markov chains have a terminal state

  • The terminal state can only transition to itself
    
  • Once the system is in the terminal state, we say that the episode has ended
    
  • Episodic processes are not of direct interest here

• Continuous Markov chains have no terminal state

  • Continue indefinitely, at least in principle
    
  • Continuous Markov chains sample probability distribution
    
  • Are ideal for estimating Bayesian posterior distributions

From Markov process to Markov chain

Markov chain comprises a number of state transitions

• Chain of \(n\) state transitions, \(\{t_1, t_2, t_3, \ldots, t_n \}\)

• Each transition has the probabilities given by the state transition matrix, \(\Pi\)

• To estimate the probabilities of being in the states we use a special case known as a stationary Markov chain

  • We will skip the technical mathematical details here
    
  • Over a large number of time steps the number of times the states are visited is proportional to the state probabilities

From Markov process to Markov chain

Start with initial state, \(\vec{x}_0\) for a continuous Markov chain:

\[\Pi\ \Pi\ \Pi\ \ldots \Pi\ \vec{x}_t = \Pi^n\ \vec{x}_t \xrightarrow[\text{$n \rightarrow \infty$}]{} \vec{p(x)}\]

• We can find the probabilities of the states without knowing the values of the transition matrix, \(\Pi\)!

• As long as we can repeatedly sample the stochastic Markov process, we can estimate the state probabilities

• This is the key to Markov Chain Monte Carlo sampling

MCMC and the Metropolis-Hastings Algorithm

The first MCMC sampling algorithm developed is the Metropolis-Hastings (M-H) algorithm; often referred to as Metropolis algorithm or the M-H algorithm.

• The M-H algorithm uses the following steps to estimate the posterior density:
  • Pick a starting point in the parameter space
    
  • Sample the posterior distribution according to the model, the product of the likelihood \(P(data|parameters)\) and prior, \(P(parameters)\)
    
  • Choose a nearby point in parameter space randomly and evaluate the posterior at this point
    
  • Use the following decision rule to accept or reject the new sample:
    • If the likelihood, \(p(data | parameters)\), of the new point is greater than the current point, accept new point
      
    • If the likelihood of the new point is less than your current point, only accept with probability according to the ratio:

\[Acceptance\ probability\ = \frac{p(data | new\ parameters)}{p(data | previous\ parameters)}\].

• If the sample is accepted, compute the posterior density at the new sample point

• Repeat sampling steps many times, until convergence

MCMC and the Metropolis-Hastings Algorithm

Eventually, the M-H algorithms converges to the posterior distribution

• M-H random sampling algorithm is far more sample efficient than naive grid sampling

• Consider that the M-H algorithm probabilistically samples the parameter space

  • Preferentially sample high density areas
    
  • Not every point on a grid

• Important properties of the Metropolis-Hastings MCMC algorithm include:

  • The M-H algorithm is guaranteed to eventually converge to the underlying distribution, but convergence can be quite slow
    
  • High serial correlation from one sample to the next in chain gives slow convergence

MCMC and the Metropolis-Hastings Algorithm

Poor convergence arrises from low sample efficiency

• Algorithm must be ‘tuned’ to ensure sample efficiency

• Tuning finds a good dispersion parameter value for the state sampling distribution

• Dispersion parameter determines the size of the jumps the algorithm makes

• Example, for Normal distribution pick the variance \(\sigma^2\)

  • \(\sigma^2\) is too small, the chain only slowly searches
    
  • \(\sigma^2\) is too big, chain has are large jumps which slows convergence

M-H algorithm example

Let’s try a simple example, find an estimate of the posterior density of a bivariate Normal distribution

\[\mu = \begin{bmatrix} .5\\ .5 \end{bmatrix}\]

\[Covariance = \begin{bmatrix} 1, .6\\ .6, 1 \end{bmatrix}\]

Random samples from a bivariate Normal distribution

M-H algorithm example

And, the marginal distributions of the variables

Marginal distributions of bivariate Normal samples

M-H algorithm example

Now, we are ready to sample these data using the M-H MCMC algorithm

The algorithm is

1. Compute the bi-variate Normal distribution likelihood
2. Initialize the chain
3. Initialize some hyperparameters statistics
4. Sample the likelihood of the data using the M-H algorithm.

MCMC samples from bivariate Normal distribution

M-H algorithm example

Plot the first 500 samples

First 500 MCMC samples from bivariate Normal distribution

Notice the long ‘tail’ on the sampled distribution from the burn-in period.

M-H algorithm example

Marginal distributions of the MCMC samples, less first 500

First 500 MCMC samples from bivariate Normal distribution

These marginals are similar to the original distribution

Convergence and sampling efficiency of MCMC

Hoe can we understand the onvergence properties of the M-H MCMC sampler

• MCMC sampling generally convergences to the underlying distribution, but can be slow

• In some pathological cases, convergence may not occur at all

• The acceptance rate and rejection rate are key convergence statistics for the M-H algorithm

  • Low acceptnce or high rejection rate are signs of poor convergence
    
  • Too few rejections, indicate that the algorithm is not exploring the parameter space sufficiently
    
  • Trade-off between these statistics is controlled by the dispersion of the sampling distribution
    
  • Unfortunately, there are few useful rules of thumb one can use

For our example these statistics are fairly good:

\[Acceptance\ rate = 0.81\] \[Rejection\ rate = 0.19\]

Evaluation of MCMC sampling

Trace plot of the samples displays the sample value of the parameter with sample number

Trace plots for two variables from the M-H algorithm

The traces show sampling around the highest density values of the parameters, indicating good convergence

Evaluation of MCMC sampling

An autocorrelation plot shows how a sample value is related to the previous samples

• The autocorrelation of the Markov chain sampling

Autocorrelation of Markov chain for the two parameters

• Notice that the autocorrelation dies off fairly quickly with lag

• Low auto correlation indicates good sampling efficiency

Evaluation of MCMC sampling

We can relate sampling efficiency to the autocorrelation of the samples

• Intuitively, uncorrelated samples provide maximum information on the distribution sampled
  • With significant autocorrelation, the new information gathered per-sample is reduced
• Effective sample size or ESS.is the ratio between the number of samples adjusted for the autocorrelation and the hypothetical number of uncorrelated samples, N

\[ESS = \frac{N}{1 + 2 \sum_k ACF(k)}\]

• ESS close to N indicates low autocorrelation and high sample efficiency

Other MCMC Sampling Algorithms

A number of other powerful MCMC sampling algorithms have been developed

• The M-H can suffer from slow convergence for several reasons
  • Generally has fairly high serial correlation and low ESS
    
  • Must ‘tune’ the state selection probability distribution
• As a result of these limitations, other MCMC sampling methods have been proposed in a quest to improve sample efficiency including:
  • Gibbs sampling
    
  • No U turn sampling (NUTS)

Gibbs sampling

Gibbs sampler is an improved MCMC sampler which speeds convergence

• Named for the 19th Century physicist Josiah Willard Gibbs; inspired by statistical mechanics
  

• Gibbs sampler samples each dimension of the parameter space sequentially in a round-robin manner
  

• M-H algorithm attempts jumps across all dimensions of the parameter space.

• Compared to the M-H, Gibbs sampling reduces serial correlation through round-robin sampling

• Update along each dimension approximately orthogonal to the preceding sampled dimension

• There are no tuning parameters since sampling is based on the marginals of the likelihood.

Gibbs sampling

The basic Gibbs sampler algorithm has the following steps:

1. For an N dimensional parameter space, \(\{ \theta_1, \theta_2, \ldots, \theta_N \}\), find a random starting point

2. In order, \(\{1, 2, 3, \ldots, N\}\), assign the next dimension to sample, starting with dimension \(1\); actual order not important

3. Sample the marginal distribution of the parameter given the observations, \(D\), and other parameter values: \(p(\theta_1|D, \theta_2, \theta_3, \ldots, \theta_N)\)

4. Repeat steps 2 and 3 until convergence

Hamiltonian MCMC

The Hamiltonian sampler was proposed in the early 1980s and uses a simple idea from classical mechanics

• Imagine that the posterior density is like a hilly landscape
  • We want to sample around the high spots, the maximum density points
• We roll an imaginary ball around the landscape to the highest potential energy points
  • These high density points attract the ball - a ‘field of attraction’
• Difficulty is that the ball will often simply roll back down the hill to a less interesting low-density regions
  • To prevent this, HMC method requires 2 hyperparameters to control stopping criteria of the ball
    
  • In practice, setting these hyperparameters proved tricky

No U-Turn Sampler

NUTS represents the state of the art in MCMC samplers

• PyMC3 package uses the NUTS MCMC algorithm.

• NUTS improves on HMC

  • The field that guides the movement of the ball is the target probability distribution
    
  • The ball is drawn towards dense (high likelihood) regions of the space
    
  • The gradient determines direction the ball rolls

• Ball is prevented from turning around

  • No U-turns!
    
  • The need for complex stopping criteria with hyperparameters is eliminated

• Why even discuss other samplers when we have the NUTS

  • While NUTS works well in many common cases, it is not guaranteed to converge
    
  • In some cases a Gibbs sampler works better

Hierarchical modeling

So far, we have only worked with models having a simple flat parameter structure

• For complex models fully employing the power of Bayesian methods, we must use hierarchical models

• Hierarchical models allow us to compute posterior distributions

  • Account for dependency in the parameters

• Build a heirarchy of parameter dependencies by factoring the joint distribution

Chain rule of probability

Use the chain rule of probability to factor a joint distribution into hierarchy

• Recall a basic relationship for conditional probability for two random variables, \(A\) and \(B\):

\[P(A,B) = P(A|B)P(B)\]

• Can expand this relationship to a joint distribution of \(n\) random variables, \(A_1, A_2, A_3, A_4 \ldots, A_n\):

\[P(A_1, A_2, A_3, A_4 \ldots, A_n) = P(A_1 | A_2, A_3, A_4, \ldots, A_n)\ P(A_2, A_3, A_4 \ldots, A_n)\]

• Continue this factorization until we reach an end point:

\[P(A_1, A_2, A_3, A_4 \ldots, A_n) =\\ P(A_1 | A_2, A_3, A_4, \ldots, A_n)\ P(A_2 | A_3, A_4 \ldots, A_n)\\ P(A_3| A_4 \ldots, A_n) \ldots P(A_n)\]

Chain rule of probability

The factorization is not unique.

• Factor the variables in any order

• For a joint distribution with \(n\) variables, there are \(n!\) unique factorizations

• Example, we can factorize the foregoing distribution as:

\[P(A_1, A_2, A_3, A_4 \ldots, A_n) =\\ P(A_n | A_{n-1}, A_{n-2}, A_{n-3}, \ldots, A_1)\ P(A_{n-1}| A_{n-2}, A_{n-3}, \ldots, A_1)\\ P(A_{n-2}| A_{n-3}, \ldots, A_1) \ldots p(A_1)\]

Factorization applied to hierarchical models

Apply the chain rule of probability to Bayes theorem

• Allows us to create hierarchical Bayesian models

• Consider an example of a multi-parameter model; estimating both the mean, \(\mu\), and standard deviation, \(\sigma\), of a Normal distribution

\[p(\mu, \sigma | D) \propto p(D| \mu, \sigma) p(\mu, \sigma)\]

• Simplify this relationship by applying the chain rule of probability for the parameter of interest, \(\mu\):

\[\begin{align} p(\mu, \sigma | D) &\propto p(D | \mu) p(\mu | \sigma) p(\sigma)\\ &\propto\ Likelihood\ *\ Conditional\ Distribution\ of\ \mu\ given\ \sigma\ *\ Prior\ of\ \sigma \end{align}\]

• Complex multi-parameter Bayesian model is transformed into a hierarchy
  • A likelihood functions, of data gien the parameter, $\(\mu\)
    
  • The conditional distribution of \(\mu\) given \(\sigma\)
    
  • The prior distibution of \(\sigma\)

The Gamma distribution

Use the Gamma distribution for the prior of the standard deviation

• Here, we parameterize the Gamma distribution by a shape parameter, \(a\) and a scale parameter, \(b\)
  • Gamma distribution is defined on the range \(0 > Gamma(a,b) > \infty\)
    
  • Desirable property for a prior, since scale parameters must have values \(>0\).

Density of the Gamma distribution with scale-2.0

Creating a MCMC hierarchical Bayes Model

We are now ready to create an hierarchical Bayes model

• It might seem some direct method of computing the posterior distribution can be easily found
  • Perhaps, but the form of the posterior is complicated
    
  • Find an approximate solution using MCMC methods
• Structure of the hierarchical model has priors for mean, \(\mu\), and standard deviation, \(\sigma\)

Hierarchical model for the posterior distribution

Creating a MCMC hierarchical Bayes Model

In mathematical terms we can define the hierarchical model as follows:

1. Prior of the scale, \(\sigma\), is \(Gamma(2, 2)\); fairly noninformative

2. Prior of the mean, \(\mu\), is \(N(0, 3)\), fairly dispersed and noninformative

3. The posterior distribution uses the priors for \(\mu\) and \(\sigma\) and a Normal likelihood

Creating a MCMC hierarchical Bayes Model

Using PyMC3 to sample the posterior distribution

def create_model(y):
    ## define a PyMC3 model object
    model = pymc3.Model()

    with model:
        # Define the prior distributions for the two model parameters, location, mu and scale
        scale = pymc3.Gamma('scale', alpha=2,beta=2)
        mu = pymc3.Normal('mu', mu=0, sigma=3.0)
        
        # The model for the posterior distribution uses the priors of the location and scale
        # to model the observations.   
        y_obs = pymc3.Normal('y_obs', mu=mu, sigma=scale, observed=y)
    return model

Evaluating the solution

Posterior distribution and trace plots for two model parameters

Evaluating the solution

PyMC3 provides quite a number of model evaluation statistics for 4 chains X 1,000 samples each

The summary shows a lot of useful information:

Summary of MCMC model performance

• Paramter means an statistics

• msce_mean is an estimate of the error arising from the MCMC sampling itself

• ESS statistics

• Gelman-Rudin statistic (Rhat) measures the ratio of the variance shrinkage between chains to the variance shrinkage within chain

Evaluating the solution

Credible intervals from MCMC sampling of model parameters

Summary

For complex Bayesian models we need a computationally efficient aproximation

• Grid sampling is inefficient

• MCMC sampling is based on Markov chains

  • Markov process is memoryless
    
  • Sampling converges to probability distribution

• Several MCMC sampling methods have been developed

  • Metropolis-Hastings (M-H) algorithm uses random sampling with acceptance probability
    
  • Gibbs sampling round robins on the dimensions of the parameter space
    
  • NUTS

• NUTS is the state of the art MCMC algorithm

  • Uses a field of attraction to the highest density parts of the parameter space
    
  • No hyperparameters to tune

Summary

Hierarchical model can be sued to represent complex relationships

• Hierarchical model based on factorization of joint distribution

  • Factorization is not unique
    
  • Choose factorization that fits the problem

• Model parameters (hyperparamters) each have a prior

• The model is sampled with MCMC algorithm

  • Find posterior of each parameter (hyperparameter)

Summary

Performance metrics of MCMC sampling

• Sample efficiency
  • Serial correlation reduces sample efficiency
    
  • ESS
• Convergence
  • Multiple chains should converge to the similar distributions
    
  • Want between chain and within chain variance to be the same